Widespread occurrence of Batrachochytrium dendrobatidis in Ontario, Canada, and predicted habitat suitability for the emerging Batrachochytrium salamandrivorans

Abstract Chytridiomycosis, caused by the fungi Batrachochytrium dendrobatidis and Batrachochytrium salamandrivorans, is associated with massive amphibian mortality events worldwide and with some species’ extinctions. Previous ecological niche models suggest that B. dendrobatidis is not well‐suited to northern, temperate climates, but these predictions have often relied on datasets in which northern latitudes are underrepresented. Recent northern detections of B. dendrobatidis suggest that these models may have underestimated the suitability of higher latitudes for this fungus. We used qPCR to test for B. dendrobatidis in 1,041 non‐invasive epithelial swab samples from 18 species of amphibians collected across 735,345 km2 in Ontario and Akimiski Island (Nunavut), Canada. We detected the pathogen in 113 samples (10.9%) from 11 species. Only one specimen exhibited potential clinical signs of disease. We used these data to produce six Species Distribution Models of B. dendrobatidis, which classified half of the study area as potential habitat for the fungus. We also tested each sample for B. salamandrivorans, an emerging pathogen that is causing alarming declines in European salamanders, but is not yet detected in North America. We did not detect B. salamandrivorans in any of the samples, providing a baseline for future surveillance. We assessed the potential risk of future introduction by comparing salamander richness to temperature‐dependent mortality, predicted by a previous exposure study. Areas with the highest species diversity and predicted mortality risk extended 60,530 km2 across southern Ontario, highlighting the potential threat B. salamandrivorans poses to northern Nearctic amphibians. Preventing initial introduction will require coordinated, transboundary regulation of trade in amphibians (including frogs that can carry and disperse B. salamandrivorans), and surveillance of the pathways of introduction (e.g., water and wildlife). Our results can inform surveillance for both pathogens and efforts to mitigate the spread of chytridiomycosis through wild populations.


| INTRODUC TI ON
Increased prevalence of fungal pathogens threatens global biodiversity, agriculture, and human health (Fisher et al., 2012;Lötters et al., 2010). The pathogenic fungi that cause chytridiomycosis (Batrachochytrium dendrobatidis (Bd) and B. salamandrivorans (Bsal)) are a well-known example of this trend. Bd affects a wide range of amphibians and is associated with population declines and even extinction of some species . The origins of Bd are enigmatic (Rosenblum et al., 2013). Southern Africa and eastern North America were initially proposed as the potential sources of this pathogen Ouellet et al., 2005), but more recent evidence suggests a common ancestral lineage in East Asia (Fisher & Garner, 2020;O'Hanlon et al., 2018). Bd now has an almost global distribution and is treated as a potentially invasive species across its contemporary range.
The related, recently discovered "salamander chytrid" Bsal is endemic to south-east Asia but recently spread outside its native range. Bsal was introduced to western Europe prior to 2004 through the pet trade (Lötters et al., 2020), and it is associated with rapid population declines of fire salamanders (Salamandra salamandra) in the Netherlands and Germany (Lötters et al., 2020;Martel et al., 2013;Spitzen-van der Sluijs et al., 2013). Populations that survived this devastating outbreak showed no subsequent resistance to Bsal (Stegen et al., 2017), suggesting that evolutionary rescue and effective prophylactic treatments (i.e., development of vaccines) are both unlikely. Further spread of Bsal is cause for concern because thousands of urodele species endemic to other regions are highly susceptible to this pathogen and would likely experience population declines if it were introduced into their ranges (Martel et al., 2014;Stegen et al., 2017). Some anuran species can also act as reservoir hosts for Bsal, potentially increasing the rate of spread following new introductions (Lötters et al., 2020;Nguyen et al., 2017).
Understanding and limiting pathogen dispersal are major objectives in the management of emerging wildlife pathogens (Biek & Real, 2010;Wilder et al., 2015). Based on the current evidence, preventing new introductions and slowing range expansion following introductions should be considered the only feasible control measures for Bsal (Stegen et al., 2017).
Species distribution models (SDMs) are a common tool used to relate environmental conditions to occurrence data for a species to predict its distribution ). In the context of pathogen management, SDMs can identify areas to prioritize for surveillance or to limit range expansion (Raffini et al., 2020). SDMs were used to predict the distribution of the "fundamental niche" of Bd in the New World and Australia (Rödder et al., 2008;Ron, 2005) and the fundamental niche of Bsal in Europe (Feldmeier et al., 2016;Lötters et al., 2020). SDMs have been used to forecast the distribution of Bd at a global scale, to quantify associated Bd-related extinction risks for amphibians under a range of climate scenarios, and to predict suitable habitat for Bsal in Mexico should an introduction occur (Basanta et al., 2019;Lötters et al., 2010;Olson et al., 2021;Rödder et al., 2009Rödder et al., , 2010. Building an SDM for an invasive species is challenging in the early stages of the invasion when data are limited (Bertolino et al., 2020;Katz & Zellmer, 2018), or when predicting suitable habitat based on occurrence data from regions with different environmental conditions. Field sampling and SDMs have been used to characterize the niche used by Bd in tropical regions (Bacigalupe et al., 2019;Catenazzi et al., 2011;Fisher et al., 2009;Kriger & Hero, 2007;Saenz et al., 2015). Models of habitat suitability for Bd in Nearctic areas have relied heavily on these data from warmer climates and have predicted low habitat suitability for Bd in northern Nearctic regions Xie et al., 2016). However, sampling in Canada detected Bd at high latitudes: 60°18′N (Schock et al., 2010); 54°73′N (D'Aoust- Messier et al., 2015); 46°81′ N (Ouellet et al., 2005). These detections could indicate a recent range expansion-but Bd was also detected in museum specimens from Quebec, Canada, collected as early as 1961 (Ouellet et al., 2005). Further sampling across Canada continues to detect Bd even with low sampling effort, and reported prevalence in samples from wild populations of Canadian amphibians ranges from 26.9% to 100% (Brunet et al., 2020;D'Aoust-Messier et al., 2015;Forzán et al., 2010;Jongsma et al., 2019;McMillan et al., 2020;Schock et al., 2010). Taken together, these studies contradict the predictions of available SDMs and suggest Bd is established and enzootic in Canada, highlighting the importance of representative sampling to inform SDMs.
Although Bsal has not yet been detected in North America, suitable climatic conditions and brisk international trade in amphibians place the continent at high risk of a Bsal invasion (Govindarajulu et al., 2017;Waddle et al., 2020;Yap et al., 2017

T A X O N O M Y C L A S S I F I C A T I O N
Community ecology; Disease ecology; Invasion ecology largely temperature-dependent (Feldmeier et al., 2016). This hypothesis was supported by experimental infection of the Nearctic eastern newt (Notophthalmus viridescens) with Bsal, which led to higher mortality at lower temperatures (Carter et al., 2021). These results shifted concern about Bsal northwards from previous risk assessments, identifying the northeastern United States and southeastern Canada as high-risk areas.
In this study, we aimed to reduce the uncertainty around the drivers and distribution of chytridiomycosis in northern climates.
We used SDMs based on local occurrence data to estimate the likely distribution of Bd in our study area (Ontario, Canada, and an island of Nunavut just off the coast of Ontario), and we used available data on suitable habitat for Bsal to identify particularly high-risk areas in the region that could be prioritized for surveillance in the case of an introduction. We conducted surveillance for Bd and Bsal in the study area across a latitudinal gradient of 15 decimal degrees that exceeds the northern range limits predicted by recent SDMs Xie et al., 2016), but within which Bd has been shown to occur (D'Aoust-Messier et al., 2015;McMillan et al., 2020). Our results are intended to inform surveillance activities and mitigation strategies to control these pathogens and their effects on amphibian diversity in temperate climates.

| Surveillance for Bd and Bsal
To obtain samples across a wide geographic area, we provided sampling kits to a broad network of government biologists and interested researchers, who opportunistically swabbed any amphibians they encountered during their work across the province of Ontario and on Akimiski Island, Nunavut (Canada). Sampling methods were approved by the Wildlife Animal Care Committee of the Government of Ontario and by Ontario Parks. Briefly, participants identified each amphibian they captured to species. They recorded potential clinical signs of chytridiomycosis, but collected swab samples from each individual regardless of clinical signs. Participants wore new nitrile gloves to handle each individual to avoid crosscontamination of samples. Samples were collected by rubbing sterile swabs (Puritan, 6" Sterile Polyester Tipped Applicator 25-806 1PD) on the ventral surface, under the forelimbs, and along the sides of the mouth. Tadpole swabbing focused on the mouth area where the fungus is able to grow on the keratinized mouthparts. Swabs were rotated while sampling to ensure all sides of the swab made contact. Samples were stored in lysis buffer, kept cool and out of direct sunlight while in the field, and then stored in fridges and/or freezers (as available) prior to laboratory analysis at the Animal Health Laboratory (AHL) in Guelph, Ontario. We used qPCR analysis of the samples to estimate distribution and prevalence of Bd and to conduct surveillance for Bsal.
The AHL facility follows strict quality control procedures, including use of a separate DNA extraction room, master mix assembly room, and amplification room. Positive control and negative control samples were included in all procedures from DNA extraction to amplification. We extracted DNA from the swab samples using a commercial DNA extraction kit (DNeasy Plant Mini Kit, Qiagen, Toronto, Ontario) and following the manufacturer's protocols. Real-time polymerase chain reaction (qPCR) assays were performed to detect Bd and Bsal, as described previously (Blooi et al., , 2016, and the qPCR assays conducted at the AHL are ISO/IEC 17025 accredited and periodically validated to ensure accuracy. To avoid cross-talk between the Cy5 and Fam fluorescent markers, we performed the duplex qPCR as two simplex qPCRs. Samples were run singly, but we confirmed inconclusive (borderline positive) samples by retesting in triplicate. We scored samples with cycle threshold (CT) values (CT) >40 as negative for the target (Bd or Bsal). This conservative approach reduces the risk of false positives, but may underestimate the true prevalence of the target pathogens. Thus, our data represent an estimate of minimum prevalence.
We tested for spatial autocorrelation among sampling locations using the join count test from the R package spdep (Bivand & Wong, 2018).

| Latitudinal distribution of Bsal detection
To explore the hypothesis of a northern range limit for Bd, we tested the association between Bd prevalence and sampling latitude using logistic regression in R 4.0.3 (R Core Team, 2020). We restricted this analysis to those species with a sample size >100 that were wellsampled across a broad latitudinal range (Lithobates sylvaticus, L. clamitans, L. pipiens, and Anaxyrus americanus).

| Species distribution models-Batrachochytrium dendrobatidis
We created six SDMs in Maxent version 3.4.1 (Phillips et al., 2006) to identify suitable habitat for Bd across our study area using our detected cases and 20 environmental variables. Maxent is an opensource program that uses presence points and a maximum entropy algorithm to estimate the probability of species occurrence when constrained by environmental conditions. Maxent has high discriminatory power (Bradie & Leung, 2017) and has been used previously to create SDMs for Bd and Bsal (Bacigalupe et al., 2019;Basanta et al., 2019;Lötters et al., 2010;Rahman, 2020;Rödder et al., 2008;Ron, 2005;Seimon et al., 2015). These studies used 19 bioclimatic variables from WorldClim 2.1 (https://world clim.org/data/index.html) (Fick & Hijmans, 2017) (Table 1) and found that many are predictors of Bd habitat suitability and chytridiomycosis severity. We also used these variables, but instead of using the 50-year averages from WorldClim, we recreated these variables for the period of our study (2014)(2015)(2016)(2017) using data from the Daily Surface Weather Data version 3 (Daymet). This NASA-sourced raster dataset includes monthly summaries of daily weather parameters across North America at a 1 km resolution (Thornton et al., 2017). We downloaded 4 years of temperature and precipitation data and used the R package dismo to recreate the 19 bioclimatic variables for each year (Hijmans et al., 2021). After averaging each variable across years, they were checked for errors using reference code available through the USGS (O'Donnel & Ignizio, 2012). To explore the effect of elevation on the distribution of Bd, we also included a Digital Elevation Model (DEM) from the Shuttle Radar Topography Mission (NASA JPL, 2013). We resampled this DEM to a 1 km grid to match the bioclimatic variables. We initially considered also including an urbanization variable, distance to roads, to account for potential human impacts on Bd occurrence. However, almost all of our samples were collected in close proximity to roads due to our opportunistic sampling design. Thus, this variable could have simply confounded the models (i.e., resulted in predictions that sites far from roads were not suitable habitat simply because our sampling was biased toward sites near roads), and we did not include this variable in our analyses.
Our six models were constructed using a combination of variable sets (Full set vs Reduced set based on a correlation analysis), and three different bias grids (a raster surface that represents sampling effort, cropped to a buffer of 0.5, 1, and 2 decimal degrees) that would restrict background point selection to address sampling bias.
A Maxent model's goodness of fit is estimated by gain (Phillips et al., 2017), which can be increased by addressing collinearity and reducing the environmental layer inputs in the model (Kramer-Schadt et al., 2013). We created a correlation matrix of variables listed in Table 1 to identify highly correlated variables (Table A1; R > .75).
We ran logistic regression among these and retained the variable that showed the strongest association with Bd presence based on point-wise values from all 1,041 sample locations. In cases where the difference was small and the choice was between variable representing averages and weather extremes, we retained the extremes as these are most relevant to niche limits Seimon et al., 2015). This process identified the nine variables used in the "Reduced" models (columns 1 and 2, potential variables into an SDM or to reduce the number of variables a priori (Feng et al., 2019). In ambiguous situations where the drivers of species and climate relationships are unknown, retaining one variable from a highly correlated cluster may cause the future range predictions to vary considerably (Braunisch et al., 2013). Therefore, we also ran full models containing all 20 variables for each of the three bias grids. Our goal was not to identify a single, best-fit model, but to evaluate the qualitative agreement among biologically informative model outputs (i.e., how robust were our results to different model parameters).
We converted all 20 Daymet variables listed in Table 1 to the same coordinate system (WGS 1984) before clipping them to a 5 km-buffered boundary of the sampling area, which included the province of Ontario and Akimiski Island (Nunavut) for processing in Maxent version 3.4.1 (Phillips et al., 2006). Maxent analyses rely on the assumption that detection probability is constant across all sites (Phillips et al., 2009), although this is frequently violated in practice (Yackulic et al., 2013). Uneven sampling effort can be addressed by thinning the number of occurrences (spatial filtering), by manipulating the background data to account oversampled areas (i.e., using a bias grid), and by reducing multicollinearity (Kramer-Schadt et al., 2013). Most of the samples were collected in southern Ontario, while the remote north-west region of the province was sparsely sampled.
We applied a spatial filter to the data by retaining only one positive For all six models, the optimal model settings were first selected using the R package ENMeval prior to running in Maxent. This package's function ENMevaluate helps fine-tune ecological niche models by running multiple iterations of user defined specifications for the regularization multiplier (RM) and Maxent feature types (e.g., linear) (Muscarella et al., 2014). The RM is a smoothing parameter that adds constraints to the model to prevent overfitting (Phillips & Dudík, 2008). We used RM increments of 0.5, between 0.5 and 2, with four combinations of linear, quadratic, product, and hinge features ("L," "LQ," "LQP," and "LQH") to create 20 candidate models using the random-fold method and 10 k-folds. We selected the models with the lowest delta Akaike's information criterion for small sample sizes (ΔAICc), a measure of goodness of fit that penalizes overparameterization (Muscarella et al., 2014). The best-supported model determined the RM and feature settings in Maxent.
We used the default ClogLog output format in Maxent, which provides an estimate of probability of presence between zero and one for each pixel. It is considered the most appropriate for estimating probability of presence, although it often estimates higher suitability values than the logistic output (Phillips et al., 2017). We partitioned Bd occurrence data into two groups (75% training, 25% testing) and tested the Maxent models with cross-validation using the default of a maximum of 10,000 randomly sampled background points. We set the maximum number of iterations to 5,000 and ran 10 replicates of each model, using the averaged suitability probabilities for our analyses.
To evaluate the individual model performance, we used receiver operating characteristic (ROC) area under the curve (AUC). A model with an AUC = 0.5 has no discrimination ability, while models with an AUC > 0.7 have adequate power of discrimination (Hosmer et al., 2013). We further evaluated the predictive accuracy of our models by calculating the true skill statistic (

| Predicting the distribution of Batrachochytrium salamandrivorans
We used two approaches to identify likely high-risk areas for the introduction of Bsal in our study area. First, we replicated the approach used by Basanta et al. (2019) in Mexico to predict likely suitable habitat for Bsal within our study area, based on the environmental conditions in areas where Bsal is native or has become established.
Second, we mapped the ranges of the species in our study area that are likely to be susceptible to Bsal, identified areas of high species richness based on those ranges, and overlaid climate variables associated with severe Bsal chytridiomycosis (Carter et al., 2021) to identify high-priority areas for surveillance in case of an introduction.
2.3.1 | Species distribution models-Batrachochytrium salamandrivorans Basanta et al. (2019) used Bsal presence data from Asia and Europe to predict where Bsal might thrive in Mexico if introduced and then compared the output to species richness in the region to identify hotspots of conservation concern. We replicated this approach for our study area using 44 Bsal occurrence records summarized in and Yuan et al. (2018). Daymet is unavailable outside North America, so we used the 19 bioclimatic variables from WorldClim 2.1 to obtain comparable environmental data for the current Bsal range and our study area (Fick & Hijmans, 2017). We selected terrestrial ecoregion boundaries from https://datab asin.org that overlapped with Bsal cases to clip the raster layers in Asia and Europe (Olson et al., 2001). We selected bioclimatic variables previously associated with Bsal prevalence (Table 2), and these eleven bioclimatic rasters were tested for multicollinearity (R > .75) as described above for the Bd models. The only difference was the use of 100 randomly sampled pseudoabsences to compare to the Bsal positives in the logistic regression tests.
To predict habitat suitability for Bsal in our study area, we constructed six SDMs based on the Asian and European occurrence data, spatially filtered so that each square kilometer only contained one occurrence. Two SDMs used all 215 cases, one with the full set of predictor variables and one using the reduced set. We also ran full and reduced models using only the occurrence data from Asia (N = 30) and only that from Europe (N = 185). We did not include elevation in these models because of the large difference in elevation between Ontario and the native range of Bsal (>3,400 m).
The R function ENMevaluate was used again to optimize the settings for each model, with RM values set to increments of 0.5 between 0.5 and 2.5, with four possible feature classes (L, LQ, LQP, LQH). We used the random-fold method for the two models with the complete dataset, and the jackknife method for the Asia-only and Europe-only models because of their smaller sample size (Muscarella et al., 2014). These outputs provided the optimized settings for the Bsal models in Maxent, which were set to 5000 iterations. A 100 km bias grid buffer around the occurrences was used to restrict the selection of 10,000 background points. Finally, we reclassified the ClogLog output of the projected suitability maps for Bsal in Ontario from high (>0.75) to no suitability (<25).

| Predicting high-risk areas for Bsal infection
We aggregated distribution data for Bsal-susceptible salamander species to quantify the distribution of species at high risk from Bsal in our were scaled from 1 to 4, and temperatures between 14°C and 22°C scaled from 4 and 1 (Carter et al., 2021). After averaging the two layers, we overlapped the result with the salamander species richness map to identify areas of high Bsal risk to native salamanders.

| Surveillance for Bd and Bsal
We collected 1,041 skin-swab samples from 18 species of frogs and salamanders across 735,345 km 2 and a latitudinal gradient of 15.07 decimal degrees (41°73′N-56°80′N), and tested these for Bd and Bsal using qPCR (Figures 1 and 2;  Figure A1).
The detected prevalence of Bd increased as the number of sampled animals increased. This trend was borderline significant using linear regression (p = .057), though the sample size was small (i.e., 18 species). Sampling locations were spatially autocorrelated (p < .0001), and most were near roads due to the opportunistic nature of the sampling.

| Latitudinal distribution
There was no linear association between Bd detection and latitude for well-sampled species (>100 occurrences/species; N = 732; p = .19), nor any apparent association for the other species (Figure 2).

| Species distribution models-Batrachochytrium dendrobatidis
We retained 86 of the 113 Bd occurrences after spatial filtering ( Figure 3a). Results from the six Bd SDMs are summarized in Table 4, and all performed best with a combination of linear, quadratic, and hinge features (LQH) with a RM of 1.5 or 2 ( Figure A2). of the study area was classified as highly suitable in the six models (Table 4). All models had omission rates below 27% (mean omission rate: 18%), and the models with the best omission rates (closest to 10%, the theoretical expectation) were the Reduced 2.0 model and Full 0.5 model, which each had omission rates of 0.11 indicating lower overfitting.
The Reduced 0.5 model was selected as the best model due a high AUC score, the highest TSS score (0.47) and the most conservative prediction of Bd habitat suitability (44.6%), although the omission rate was moderate (0.18). Figure 3a,b present the mapping output for the Reduced 0.5 model and the area classified as suitable using the 10th percentile threshold, respectively. Maps based on the other five models are available in the Appendix A ( Figure A3, Figure A4). Results were similar, indicating that these models were qualitatively robust to variations in model inputs and settings.
Three variables contributed to 68.5% of the variability and drove performance of the Reduced 0.5 model (Table A3) gesting that BIO18 improves model fit ( Figure A6).

| Species distribution models -Batrachochytrium salamandrivorans
Our Bsal SDMs all predicted low habitat suitability across our study area, based on extrapolation from Bsal occurrences in Asia and Europe. They predicted very little suitable habitat in southern Ontario, and the highest suitability fell along the coast of Hudson's Bay, at the northern limit of the province. This result was clearly not biologically meaningful and is most likely caused by our attempt to predict future distributions based on data from regions with very different climatic conditions. We therefore discarded our predictive Bsal SDMs and did not pursue this approach further, although the results are presented in the Appendix A ( Figure A7) for transparency. Each swab was also tested for Bsal, but that pathogen was not detected in any sample. a Five of these individuals were identified as possible laterale-jeffersonianum hybrids.

| Predicted high-risk areas for
were assigned to an area on the north shore of Lake Superior, followed by all of southern Ontario, and part of the coast of James Bay ( Figure 4b). Over 90,123 km 2 , or 9.8% of the province, had scores ≥2.5 indicating moderate to high vulnerability (Figure 4c). By overlaying areas with ≥4 vulnerable salamander species and areas with temperature-derived risk scores ≥2.5, we identified high-risk areas covering 60,530 km 2 (6.6% of Ontario) that could be prioritized for Bsal surveillance or mitigation if an introduction occurs (Figure 4c).
The two layers overlapped substantially; 50.4% of the area containing ≥4 vulnerable species also fell within the higher temperature vulnerability zone. F I G U R E 4 (a) Salamander richness map of five Urodele species likely susceptible to Bsal. All point data were collected between 1950 and 2020 and are a combination of our records and GBIF occurrences (GBIF, 2020). (b) Temperature-based Bsal risk scores following Carter et al. (2021). (c) Area of overlap between high species richness (>4 species) and higher temperature risk scores (>=2.5), which covers 60,530 km 2

| DISCUSS ION
area (See also D'Aoust- Messier et al., 2015;Fisher & Garner, 2020;Ouellet et al., 2005). We were encouraged that we did not detect Bsal in our sampling. Although our Bsal SDMs proved uninformative, they nicely illustrate the risks (as outlined in the Introduction) of constructing predictive SDMs that rely on occurrence records from very different climatic regions. In lieu of SDMs, we were able to more coarsely identify high-risk areas for Bsal invasion by overlaying temperature-derived risk scores with range maps of vulnerable species. These results can inform Bsal surveillance in Ontario, and once Bsal occurrence data from more representative locations become available, it will be possible to generate more informative, predictive SDMs. Ontario, and Atlantic Canada) detecting Bd prevalence ranging from 13.1-100%, (Forzán et al., 2010;Longcore et al., 2007;McMillan et al., 2020;Ouellet et al., 2005). In our study, the estimated prevalence of Bd should be interpreted as estimates of minimum prevalence. Surveillance based on qPCR can only detect fungal loads above the detection threshold, missing cases where pathogen load is low. We also used skin swabs that are non-invasive and easy to replicate, but are less sensitive for detecting Batrachochytrium than tissue samples (Schock et al., 2010). These factors all lead to underestimates of true prevalence. Likewise, the results of our SDMs represent minimum estimates of suitable habitat for Bd, because we used a relatively small number of presence points (N = 86; Proosdij et al., 2016). Some areas were also undersampled as they do not have road access and are more difficult to reach. Nondetection of Bd from areas such as the Hudson Bay coast probably reflect low power to detect the pathogen with small sample sizes, rather than lower prevalence in the north. Despite these limitations, we detected Bd as far north as Akimiski Island, which had an average temperature of −1.25°C over our study period, and our most conservative SDM for Bd predicted that suitable habitat extends into the sub-arctic (54°38′N). Bd grows slowly or not at all at temperatures <4°C, but can overwinter at lower temperatures (Piotrowski et al., 2004).
Zoospores have slower rates of growth at lower temperatures, but exhibit longer activity, which may increase encounter rates with amphibian hosts (Voyles et al., 2012;Woodhams et al., 2008). Thus, the thermal sensitivity of Bd allows it to tolerate conditions common to our study area, contrary to previous niche models that identified the boreal and northern parts of North America as unsuitable habitat Xie et al., 2016). Including regionally relevant occurrence data in our Bd SDMs resulted in a much broader predicted niche and highlighted the ubiquity of Bd in our sampling area.
Ironically, lack of regionally relevant occurrence data also prevented us from predicting suitable habitat for Bsal. We were still able to broadly identify areas where salamanders would be most at risk from a potential Bsal introduction. However, the maps we generated through this exercise represented a tool to highlight priority areas for surveillance in the event of a Bsal introduction, rather than a prediction of the full, available niche for Bsal in Ontario. Use of these maps to plan surveillance for Bsal must consider three important caveats. First, the species richness layers were created using broadscale range maps that did not capture local variation in the distribution of each species. Second, our risk scores for Bsal were derived from experimental infection of N. viridescens (Carter et al., 2021).
Use of these risk scores required an explicit, untested assumption:  (Forzán et al., 2010). Periods of high Bd prevalence in sampled roadkill were also not correlated with documented die-offs in Maine, USA (Longcore et al., 2007). Chytridiomycosis severity varies with host susceptibility and strain virulence (Eskew et al., 2018;Gahl et al., 2012), and the strain or strains present in our study area may simply be less virulent than those occurring elsewhere. Continued surveillance can assess the potential emergence of more virulent strains, such as that observed in a population of Cascades frogs (Rana cascadae) in northern California (Piovia-Scott et al., 2015). Nevertheless, although Bd has had severe impacts on many amphibian populations, amphibians in our study area appear to currently co-exist alongside it with only intermittent outbreaks (Kilpatrick et al., 2010).
In contrast, the impacts of a Bsal introduction on salamander in our study area could be severe. Bsal grows optimally at 10-15°C, can survive between 5 and 25°C, and has a lower thermal preference than Bd . We identified southern Ontario and two regions farther north as likely high-risk areas in the event of a Bsal introduction. Southern Ontario is one of the most densely populated regions in Canada, and wildlife in this area are already threatened by habitat loss and fragmentation (Paterson et al., 2021).
The large cities in southern Ontario represent potential trade routes for introducing Bsal into the country (Stephen et al., 2015), highlighting international trade restrictions as an important policy tool (Auliya et al., 2016;Stegen et al., 2017). As of May 31, 2017, importation of salamanders into Canada was prohibited without a permit (Government of Canada, 2017). However, anurans can carry Bsal without exhibiting clinical signs of disease  and can transmit the pathogen to salamanders with fatal results (Stegen et al., 2017). Expanding trade restrictions to include anurans could further protect native amphibians.
In the era of broad-scale ecological problems that require in-  (Stephen et al., 2015), and their recommendations remain relevant. These include not only amending import policies but also applying obligatory testing for Bsal, early detection and containment, and increasing awareness of the pathogen within the pet trade community. The current ban on salamander importation is encouraging, but the potential importation of Bsal on anuran hosts remains a concern. Policy reforms, continued surveillance, and international cooperation will be required to prevent potential introductions into North America and protect native amphibians. We also recognize that present-day Ontario and Akimiski Island exist on the traditional territories of multiple Indigenous nations.

ACK N OWLED G M ENTS
Conservation of at-risk amphibians in our study area is made possible, in part, by the historical and ongoing land stewardship of Indigenous communities that support the persistence of healthy wildlife populations.

This study was funded by the Government of Ontario and the
Canadian Wildlife Health Cooperative. The authors declare no conflicts of interest.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data Badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at https://doi.org/10.17605/ OSF.IO/7XNRE.   These species are confirmed or considered potentially susceptible to Bsal (Martel et al., 2014).

TA B L E A 2
Species in dataset obtained from our occurrence records and the Global Biodiversity Information Facility (GBIF, 2020)

F I G U R E A 1 Deceased
Lithobates pipiens that tested positive for Batrachochytrium dendrobatidis near Peterborough, Ontario, Canada

F I G U R E A 2
Results of ENMevaluate in R for selecting Regularization Multiplier (RM) for the (a) Reduced Bd Models and (b) Full Bd models will 20 variables. Lower ∆AICc values indicated a better model F I G U R E A 3 (a) Reduced Model outputs of predicted habitat suitability of Bd using the ClogLog distribution using three different background extents for the bias grid (0.5, 1, and 2 decimal degrees). (b) The 10th percentile training presence ClogLog threshold for each background bias grid. Shading indicates the area of habitat suitability. The reduced models were created with nine predictor variables listed in Table 1 (a)

(b)
F I G U R E A 4 (a) Full Model outputs of predicted habitat suitability of Bd using the ClogLog distribution using three different background extents for the bias grid (0.5, 1, and 2 decimal degrees). (b) The 10th percentile training presence ClogLog threshold for each background bias grid. Shading indicates the area of habitat suitability. The full models were created with 21 variables listed in

F I G U R E A 7
Reduced Model output of forecasted habitat suitability for Bsal using the ClogLog distribution using all available occurrence data for Bsal from (a) Asia and Europe, (b) Asia only, (c) Europe only. (d) Classification of Bsal habitat suitability projection based on the model using only Asia and all 29 variables. (<0.25 = no suitability, 0.25-0.5 = low suitability, 0.5-0.75 = moderate suitability, >0.75 = high suitability), based on the methods of Basanta et al. (2019). These results are unlikely to be accurate or informative for making meaningful predictions about Bsal habitat suitability in our study area. We therefore excluded these results from further analyses or interpretation in the main text. However, we have chosen to include them in this Appendix because they highlight the challenges in extrapolating Species Distribution Models among regions with vastly different climates and conditions